!> Only structured meshes supported so far
module mesh_mod

use global_mod, only: ink, rlk

use point_mod,  only: PointObj

use memory_mod, only: Allocate

implicit none

real(kind=rlk)    :: xLimits(2) = (/0.0,0.0/)
real(kind=rlk)    :: yLimits(2) = (/0.0,0.0/)
real(kind=rlk)    :: zLimits(2) = (/0.0,0.0/)
integer(kind=ink) :: xCells = 0
integer(kind=ink) :: yCells = 0
integer(kind=ink) :: zCells = 0

namelist /main_mesh/ xLimits, yLimits, zLimits, xCells, yCells, zCells

!> Mesh object description
type MeshObj
  logical :: zinit = .true.
  integer(kind=ink) :: nNodes = 0
  integer(kind=ink) :: nCells = 0
  integer(kind=ink) :: nCells_x = 0
  integer(kind=ink) :: nCells_y = 0
  integer(kind=ink) :: nCells_z = 0
  integer(kind=ink) :: nNodes_x = 0
  integer(kind=ink) :: nNodes_y = 0
  integer(kind=ink) :: nNodes_z = 0
  real(kind=rlk), pointer :: axis_x(:)  => NULL()
  real(kind=rlk), pointer :: axis_y(:)  => NULL()
  real(kind=rlk), pointer :: axis_z(:)  => NULL()
  type(PointObj), pointer :: pt_list(:) => NULL()
end type

!> Mesh creation
interface Mesh
    module procedure Mesh_CreateFromDeck
end interface

contains

!> Create a mesh from scratch
type(MeshObj) function Mesh_CreateFromDeck() result(mesh)
  mesh%nCells_x = xCells
  mesh%nCells_y = yCells
  mesh%nCells_z = zCells
  mesh%nCells   = max(1,mesh%nCells_x)* &
                  max(1,mesh%nCells_y)* &
                  max(1,mesh%nCells_z)
  mesh%nNodes_x = mesh%nCells_x+1
  mesh%nNodes_y = mesh%nCells_y+1
  mesh%nNodes_z = mesh%nCells_z+1
  mesh%nCells   = mesh%nNodes_x*mesh%nNodes_y*mesh%nNodes_z

  call Allocate(mesh%axis_x, mesh%nNodes_x)
  call Allocate(mesh%axis_y, mesh%nNodes_y)
  call Allocate(mesh%axis_z, mesh%nNodes_z)

  !write(6,*) "BUILDING MESH AXES"
  if (xCells>0) call Axis_Build(mesh%axis_x, xLimits)
  if (yCells>0) call Axis_Build(mesh%axis_y, yLimits)
  if (zCells>0) call Axis_Build(mesh%axis_z, zLimits)
end function Mesh_CreateFromDeck


subroutine Axis_Build(xx, xrange)
  real(kind=rlk), intent(inout) :: xx(:)
  real(kind=rlk), intent(in)    :: xrange(2)
  integer(kind=ink) :: ii, nx
  real(kind=rlk)    :: dx
  nx = size(xx)-1
  dx = (xrange(2)-xrange(1))/nx
  xx(1) = xrange(1)
  do ii=1,nx; xx(ii+1) = xx(ii) + dx; enddo
end subroutine Axis_Build

end module mesh_mod

